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SUMMARY 


A geometrical nonlinear analysis “based on an alternative definition of 
strain is presented. Expressions for strain are obtained by con5)uting the 
change in length of the base vectors in the curvilinear element coordinate sys- 
tem. The isoparametric element formulation is assumed in the global Cartesian 
coordinate system. The approach is based on the minimization of the strain en- 
ergy, and the resulting nonlinear equations are solved by the modified Newton 
method. Integration of the first and second variation of the strain energy is 
performed numerically in the case of two- and three-dimensional elements . Ap- 
plication is made to a simple long cantilever beam. 


INTRODUCTION 


The nonlinear finite element formulation described here represents a part 
of the development of the BASIS finite element analysis system (ref. l). The 
basic idea was to combine linear and nonlinear behaviour in order to deal more 
efficiently with structural analysis. Thus, nonlinear elements had to be devel- 
oped which fit into an existing system without loss of general validity. The 
Lagrange formulation has therefore been chosen (ref. 2), expressing the dis- 
placement variables directly in the global Cartesian coordinate system, 
although locally, for practical purposes, a skew coordinate system may be pre- 
scribed. The advantage of this formulation is the numerical method arising from 
it. By adopting the modified Newton method (ref. 3 » a clear solution process 
has been chosen, therefore minimizing errors due to a wrong understanding of 
nonlinear behaviour. The strains can be adapted to the nonlinearity of the 
problem. This feature may considerably reduce the computational effort. 


* Research sponsored in part by the Swiss National Science Foundation (contract 
No. 2267-079) and the National Swedish Board for Energy Source Development 
xinder contract No. 5061-OI2. 


251 



A set of one-, two- and three-dimensional elements (in terms of curvilin- 
ear coordinates) can be derived from this formulation (Fig. 1). The elements 
presented here are based on the isoparametric approximation (ref. 3). One-di- 
mensional elements are well suited to deal with cable and truss problems 
(ref. 5)* Membrane and volume elements are mainly used to model elastic com- 
posite structures. 


GOVERNING EQUATIONS 


Consider a point P of the undeformed body (Fig. 2). The corresponding 
cartesian coordinate vector {r} may be expressed as a function of the body’s 
curvilinear coordinates 0 

n 


{r} = I {r}.>F.(0“) 
i=1 


1 1 


( 1 ) 


where are the interpolation functions of the corresponding nodes. The base 
vectors at P are obtained by deriving {r} with respect to the coordinates 0 


9{r> 
a ^ 90' 




(2) 


For practical purposes, Eq. 2 can be rewritten 

{g}^ - (3) 

[A]ot contains all functions and {s} all Cartesian components of the node 

vectors {r}. . 


In order to define the strain energy the metric tensor at P 

° 0 = 1.3; B = 1,3 

and the infinitesimal volume 

dV = det(g^g) ‘d0^ • d0^ • d0^ 

are needed. In the deformed state the base vectors become (Fig. l) 

^ n 


Eq,. 6 may be rewritten 


{G} = ^ I ({r}.'!'. +{u}.¥.) 

a 90 . . 11 11 

1=1 


{G} = [A] * {d} + {g} 

a a a 


(h) 

(5) 

( 6 ) 
(7) 


where {d} contains, similar to {s}, all Cartesian components of the node 
displacements. The metric tensor of the deformed body at P becomes 




( 8 ) 
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( 9 ) 


Substituting Eqs. 3 and 7 in Eq. 8, 


with 




[Dl^g = W; • U]^ 


( 10 ) 


The strains are now defined as the relative change of the base vectors ; 


or 


'a3 


'a3 




G + G „ + G_ + Gqq 
gg a3 3a 33 

l_ ®gg ^ ^®g3 ®36 J 


- 1 


Eq. 12 can be expressed in terms of {d} 


'g3 


' {d}'^[D]^g({d} -i-2{ s}) 

®g3 


I 

■ 2 


' + 1 


- 1 


with 


[D] Q = [D] + [D] - + [D]„ + [D]q_ 

ap Ota a3 pot 33 


«a6 ‘ %c =^8a6 * %B 


( 11 ) 


(12) 


(13) 

(li;) 

(15) 


Eq. 13 defines the strains e ^ used to derive the numerical equations from the 

r • 0^6 

strain energy. 

Note here that if the deformations are limited in size the root in Eq. 13 can 
be expressed as a series expansion according to 




(16) 


Retaining only the first two terms of Eq. l6 leads to the quadratic approxima- 
tion _ 

^ {dK[D] _({d}+2{s>) 

^ <’T> 

from which the linear solution may be obtained. 

Assuming linear elastic behaviour, the variation of strain energy leads directly 
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to the nonlinear equilibrium equations of the element. 

For ease of formulation the strains are rewritten in vector form, {e} = 

9^22 9 ^3 3 9 ^12 9^13 9 ^ 23 }^ Furthermore, as the stress-strain relations are 
mainly formulated in a rectangular coordinate system the strain is transformed 

{y} = [T]{e} (18) 


The transformation matrix [T] is constant if the scalar products of the base 
vectors in the undefonaed and in the deformed system are the same. 


The strain energy density then becomes 



dU = ^ ({a}^-{ao}^)[E]({Y> -{Yo>)<iV 

(19) 

or 

dU = -^ ({e}^[T]^[E] -{cro}^)[T]{e}dV 

( 20 ) 


Eq. 20 is then integrated over the element. For one-dimensional elements an 
analytical solution is possible (ref. 5)5 whereas for two- and three-dimen- 
sional elements numerical integration must be performed. Using Gaussian inte- 
gration the strain energy becomes 


m 


1=1 


( 21 ) 




where w designates the weight factor at the corresponding point 0 . 

X K, 

The set of • nonlinear equations for one element is directly obtained from the 
first variation of the total potential energy with respect to the global vari- 
ables 


9 U 


q. (u. ) = 
1 J 


= f. 


( 22 ) 


{f} denotes the vector of external loads. It is worth mentioning here that in 
the case of dynamic analysis the equilibrium equations can also be obtained by 
applying Hamilton’s principle and a suitable operator for discretization in 
time . 


After rearrangement, Eq. 22 becomes 



{q} = I ({£}J*[E] 

k =1 ^ 


(23) 

or 

m 

(a) = I 

k=l 


( 24 ) 

with 

[P]'' = 

a = 1 , 3 ; 3 = 0 ^ 9 3 

(25) 
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and 


(26) 


= [D]„g-(W+{s}) 

®a6 

In order to solve Eq. 22 “by the modified Newton method the Jacobian of {q} is 
needed : 

= [^] 

^ 1 J 


®aB 

For a set of several elements the stiffness matrices [K] are assembled to the 

global stiffness matrix and factored. Using the modified Newton method the n-th 

iteration becomes _ 

{d}^-{d} = 6-[K]/({f} -{q} ) (28) 

n+1 n Z n 

[K] being factored at iteration step Z ^ n. The relaxation factor is computed 
using an extrapolation method for each iteration (ref. 5). Thus the final load 
level is reached stepwise, and for each new load step the solution is extra- 
polated quadratically . Convergence criteria are based on the Euclidean and 
maximum norms ( ref . 6 ) . 

The computational procedure is essentially the same for all element types . How- 
ever, for one-dimensional elements {q} and [K] can be determined analytically 
in terms of u. , and the transformation [T] is not necessary. The [T] matrices 
for two- and 4hree-dimensional elements are listed in the Appendix. 

During the first assembly of [K] the [D] , [T] , [E] matrices as well as the 
geometry parameters g^^g and g are coii5)uted once for each Gaussian point and 
reused for further con 5 )utations of [K] and {q} . For each iteration step {e} 
and [P] are evaluated at each Gaussian point and the global {q} vector as- 
sembled. The global load vector contains not only external forces but also the 
first variation coming from linear elements . In fact , their contribution has to 
be evaluated only once for each load step, thus reducing the computing time. 


a=1,3;6=a,3 


(27) 


m 


= I 


k=1 
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NUMERICAL RESULTS 


An application of the theory is demonstrated using results obtained by the 
BASIS computer program (ref. l). A long thin cantilever beam (Fig. 3) has been 
idealized by 8-noded membrane elements. Its length is 1000 mm, height 10 mm and 
thickness 1 mm. The elastic modulus is 3000 N/mm^ and the Poisson ratio 0.36. 
The beam is subjected to a load case consisting of a variable load Q at the 
node A and a case with constant Q and variable compression load P. Since the 
strains remain smaller than one in this example differences in the results are 
not detectable when using Eq. 17 instead of Eq. 13. The convergence criterion 
during iteration (Eq. 28) is based on the relative change of the displacement 
vector norm and has been set to 0.0001. The active load is applied stepwise. If 
convergence is rapid, the step is automatically increased. However, for highly 
nonlinear problems, it is preferable to recompute and refactor the stiffness 
matrix when the iteration diverges rather than to decrease the load step. The 
first and second variation of U has been computed using 2 by 2 Gaussian inte- 
gration. Fig. h shows the load-displacement function at node A for transversal 
loading, and the same function for the divergence problem is exhibited in Fig.5* 
Note the good behaviour of the two-element approximation even for large non- 
linearities . 

To conclude, it should be mentioned that the problems currently being in- 
vestigated include the influence of the integration order on numerical accuracy 
and convergence behaviour as well as the nonlinear creep of structures . 


SYMBOLS 


Vectors are symbolized by {}-brackets and matrices by []-brackets. []^ 
means transposed matrix, inverted or factored matrix, and ^ stands for 

dyadic product. Greek indices refer to the curvilinear coordinate system. 

ot 

[A]^ matrix of form functions of curvilinear coordinate 0 

[D]^p product of form function matrices (Eq. 11) 

{d} nodal displacement vector (global Cartesian components) 


{e} 


strain vector. 


contains components e ^ 

Otp 


[E] elasticity matrix 




a6 

{g} 


a 


metric tensor of the undeformed body 
metric tensor of the deformed body 

ot 

base vector of the iinde formed body, coordinate 0 
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g 

{y} 

{yo> 

[K] 

[p] 

1 

{q} 

{r}. 

1 

{s} 

{a} 

{ao> 

[T] 

0“ 

{u}. 

1 

u 

X 

^a3 


base vector of the deformed body 

determinant of the metric tensor 

strain vector in the local Cartesian system 

initial strain vector in the local Cartesian system 

second variation of the strain energy (stiffness matrix) 

matrix defined in Eq. 26 

interpolation f mictions for geometry and displacements 

first variation of the potential energy 

global Cartesian coordinate vector of node i 

vector containing all Cartesian node coordinates 

stress vector in local Cartesian system 

initial stress vector in local Cartesian system 

transformation matrix relating {e} to {y} (see Appendix) 

curvilinear element coordinates 

global Cartesian displacement at node i 

strain energy 

weighting factor for Gaussian integration 
orthogonal Cartesian coordinates 
expression defined in Eq. 13. 
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APPENDIX 


The strain is alvays transformed to a Cartesian coordinate system defined 
by {v} = {r }2 ”{r}i for tvo-dimensional elements and by {v} and {r }4 for three- 
dimensional elements. For membrane elements, [T] becomes 


[T] 


1 

2 


{ 1 + cos2d) 

a 


1 - cos2d) 

a 


2sind) cosd) } 
a a 


a = 1,3 


{v} • {g} 


a 


COS(f) = 

“ |{v}| • |{g} 


a' 


For volume elements the six strain components are transformed similarly. 
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